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The  Software  Scene  in  The  Extraction  of  Eigenvalues 
from  Sparse  Matrices 

B.  N.  Parlett 


ABSTRACT 

The  class  of  users  of  programs  In  this  area  is  diseased  and 
split  into  the  intensive  and  the  sporadic  subsets.  Available 
software  for  each  group  is  reviewed  and  some  current  develop¬ 
ments  are  outlined. 

This  essay  arose  from  a  talk  delivered  at  the  Sparse  Matrix 
Symposium  held  at  Fairfield  Glade,  Tennessee  in  October.  1982. 
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I.  Introduction 


Numerical  analysts  see  their  task  as  the  development  and  analysis  of 
methods  for  computing  quantities  which  arise  in  science  and  engineering.  Some 
assume  that  any  good  algorithm  will  soon  find  a  user.  Who  knows  if  that  assump¬ 
tion  is  warranted?  In  any  case  this  essay  will  begin  by  taking  a  step  back  from 
software  to  consider  who  wants  eigenvalues  and  eigenvectors  of  sparse  matrices 
of  large  order.  It  is  encouraging  to  see  that  some  computing  sendees  are 
automatically  monitoring  the  use  of  programs  in  their  libraries-  (e.g.  S&ndia  Labs 
in  Albuquerque).  Such  practices  will  help  at  least  a  few  people  assess  the  actual 
usage  of  all  those  carefully  written  subroutines. 

Our  investigations  though  limited,  all  point  to  an  important  distinction 
between  users  that  goes  a  long  way  towards  explaining  the  present  state  of 
affairs  in  software  development.  That  is  the  gist  of  Sections  3  and  4. 

This  essay  may  well  be  read  by  those  who  are  not  specialists  in  software  for 
eigenvalue  calculations  and  so  there  is  a  section  which  .tells  the  EISPACK  story. 

The  very  existence  of  this  set  of  subroutines  is  intriguing  and  it  continues  to 

» 

influence  the  development  of  mathematical  software.  Our  account  is  too  brief  to 
do  justice  to  its  subject  but  it  provides  essential  background  material. 

N 

In  Sections  6,  7,  8  we  discuss  the  programs  which  were  readily  available  in 
September  1982  and  then,  in  Section  9  we  describe  some  programs  that  are  still 
under  development. 

Terminology 

It  is  vital  to  science  that  one  not  be  more  precise  tjban  is  necessary  for  the 
purpose  in  hand.  In  this  essay  words  like  small,  large,  and  sparse  play  a  key 
role.  Yet  their  meaning  will  be  relative  to  the  user’s  computing  system. 

We  say  that  an  n  x  7k  matrix  is  small  if  two  n  x  n  arrays  can  be  held  in  the 
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fast  store,  in  addition  to  any  programs  that  are  needed,  otherwise  the  matrix  is 
large.  This  binary  distinction  should  not  be  pressed  too  hard  or  it  will  crumble. 
The  purists  say  that  a  matrix  is  sparse  if  it  has  only  a  few  nonzero  elements  in 
each  row,  not  more  than  50  say.  On  this  definition  a  matrix  with  90%  of  its  ele¬ 
ments  zero  would  not  be  sparse  because  (1/  10)n2  >  50n  for  large  enough 
n(n  >  500) . 

Our  definition  suggests  that  it  is  reasonable  to  use  similarity  transforma¬ 
tions  to  help  find  the  eigenvalues  of  small  matrices.  Our  definition  also  suggests 
that  one  should  cherish  the  zero  elements  of  sparse  matrices.  Where  does  that 
leave  large,  nearly  full  matrices?  Answer  with  the  chemists  -  see  Section  4. 

The  exploitation  of  sparsity  in  the  execution  of  triangular  factorization  has 
given  rise  to  a  valuable  and  vigorous  research  area  called  sparse  matrix  technol¬ 
ogy.  The  reriews  of  Duff  describe  it  well.  [Dull.  1982].  It  turns  out  that  this 
technology  is  not  directly  relevant  to  eigenvalue  problems.  There  are  two 
classes  of  methods:  those  that  employ  explicit  nondiagonal  similarity  transfor- 

i 

mations  and  those  that  do  not.  That  is  all.  If  it  is  attractive  to  use  explicit  simi- 

v 

larity  transformations  then  please  yse  them.  Our  concern  here  is  with  the  other 
class,  whatever  the  character  of  the  matrix 

This  is  the  place  to  mention  that  G.W.  Stewart  gave  a  concise  review  of 
numerical  methods  for  sparse  eigenvalue  problems  in  [Duff  Jk  Stewart,  1976]. 
That  covers  the  background  of  most  of  the  software  we  discuss  here.  Since  then 
there  have  beerf  important  refinements  to  the  Lanczos  algorithm.  Both 
Davidson'*  method  and  the  ideas  sketched  in  the  section  on  recent  develop¬ 
ments  are  of  more  recent  vintage. 

2.  EISPACK 

EISPACK  is  an  organized  collection  of  some  40  FORTRAN  subroutines  which 


compute  some  or  all  eigenvalues  of  a  given  matrix  with  or  without  the 
corresponding  eigenvectors.  This  large  number  of  subroutines  reflects  the 
yearning  for  efficiency.  E1SPACK  has  special  techniques  for  real  matrices,  for 
symmetric  matrices,  for  tridiagonal  matrices  (ay  =0  if  k  -jfl>  1)  .  and  for 
upper  Hessenberg  matrices  (ay  =  0  if  i  —  j  >  1)  . 

The  first  issue  of  EISPACK  appeared  in  1974.  The  package  was  distributed 
by  the  code  center  of  the  Argonne  National  Laboratory,  Argcnne.  Illinois,  and  the 
indispensable  EISPACK  GUIDE  was  published  by  Springer-Vsrlag,  New  York.  To 
complete  most  eigenvalue  calculations  it  is  necessary  to  invoke  more  than  one 
subroutine.  For  example,  a  nonsymmetric  matrix  B  might  be  first  "balanced" 
by  a  diagonal  similarity  transformation,  B  ■*  DBD~ 1  =  :C  ;  next  C  is  reduced  to 
Hessenberg  form.  C  -*  P*  CP  =  :H  with  P  orthogonal  and  then  H  is  reduced 
to  triangular  form  T  (to  within  working  accuracy),  H  -*  Ql  HQ  -  :T  with  Q 
orthogonal.  The  eigenvalues  of  B  are  found  on  the  diagonal  of  T . 

Such  detailed  knowledge  of  the  process  was  more  than  some  users  wanted 
to  learn.  To  meet  this  criticism  some  members  of  the  EISPACK  team  wrote  an 
interface  called  EISPAC  which  allowed  the  user  to  specify  the  type  of  matrix  and 
the  quantities  to  be  computed.  The  interface  then  chose  the  appropriate  sub¬ 
routines  and  called  them  with  the  right  parameter  values.  Unfortunately  EISPAC 
was  only  available  on  IBM  systems. 

A  second  edition  of  EISPACK  appeared  in  1977  and  a  third  is  scheduled  for 
early  1963.  Each  edition  removed  blemishes  found  in  some  subroutines  and 
added  new  programs.  We  enlarge  on  this  topic  further  on. 

What  many  users  of  EISPACK  do  not  know  is  that  it  represented  a  coopera¬ 
tive  effort  by  most  of  the  world's  experts  in  matrix  computations.  It  seems  safe 
to  assert  that  not  (me  of  the  people  involved  in  EISPACK,  or  its  antecedents, 
became  rich  in  the  process.  What  a  contrast  to  the  software  produced  for 
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structural  analysis!  There  the  temptation  to  form  a  company  and  sell  the 
■oft wore  for  profit  was  great.  Consequently  there  are  enumerable  rival  pack¬ 
ages  which  preside  similar  services.  The  best  known  are  NASTRAN.  STRLDL 
SAP.  MARC.  ADOUL  Tins  situation  provides  a  race,  because  m.  enable  of  com- 
petition  not  leading  to  superior  performance.  Vas  the  cooperation  and  sharing 
which  characterised  the  numerical  analysts  doe  to  their  admirable  sense  a f 
values  or  to  the  absence  of  a  strong  enoqgfa  industrial  demand  for  pore  eigen¬ 
value  codes?  The  solution  of  this  conundrum  is  left  to  the  reader. 

Another  remarkable  aspect  of  QSPACK  is  that  its  first  version  defiberaleiy 
eschewed  the  production  of  new  algorithms.  The  goal  was  “simply'"  to  translate 
into  FORTRAN  same  of  the  ALGOL  programs  in  the  famous  “Handbook  for 
Automatic  Computation,  wnhana  2.  !««-»■  Algebra**,  by  J.  H.  WHkinson  and  C. 
Rrisr.h  The  authors  ware  tba  acknowiedgwd  leaders  in  the  art  of  matrix  compu¬ 
tations  The  Handbook  appaared  in  1971  and  represented  the  tndy  international 
cooperation  mentioned  earlier.  Not  all  the  programs  were  written  by  Wilkinson 
or  Rsmsch  but  each  contribution  was  carefuky  scrutinised  by  them  and.  more 
often  than  not.  improvements  were  incorporated  into  the  final  wrawu  More¬ 
over  moat  of  the  programs  had  already  been  ptddisbed  in  Numariacfae  Mathema- 
tik  and  so  had  been  refereed,  and  tasted,  and  available  in  the  public  domain 

Why  mention  these  details?  Because  the  exercise  of  translating  tzese 
ALGOL  programs  into  PCRTRAX  n«  far  from  trivial  ana.  incredible  as  .*.  may 

seem,  blemishes  were  found  in  a  number  of  the  handbook’s  program*.  Of 
course,  the  QSPACK  team  was  aiming  high.  It  was  not  enough  to  bavw  portable 
subroutines  (far  the  1906  stmsdmd  version  of  FORTRAN)  but  these  programs 
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not  present,  can  have  a  strong  effect  on  the  time  needed  to  execute  a  matrix 
computation.  Fortunately  it  was  possible  to  express  the  algorithms  so  that  they 
would  perform  satisfactorily  on  all  the  major  systems.  This  is  the  nature  of  work 
in  mathematical  software. 

The  next  aspect  of  EISPACK  that  should  be  remembered  is  the  massive 
effort  at  testing  the  programs,  in  an  adversary  manner,  before  their  release. 
About  IS  computing  groups  in  the  U.S.A.  and  Canada  agreed  to  install  and  test 
the  trial  programs  they  received  from  Argonne.  As  bugs  and  blemishes  were 
uncovered  there  were  several  iterations  on  the  programs. 

A  nasty  thought  must  be  stirring  in  the  reader's  mind.  How  can  an  algo¬ 
rithm  be  published  in  Numerische  Mathematik  after  careful  refereeing,  be  sub¬ 
ject  to  scrutiny  and  testing  by  Wilkinson  and  Reinsch  prior  to  inclusion  in  the 
Handbook,  be  translated  into  FORTRAN  by  the  EISPACK  team  with  close  attention 
to  detail,  be  tested  by  various  sites  charged  with  that  duty,  and  yet  retain  a  bug 
or  even  a  blemish?  Either  the  people  involved  in  this  development  are  not  truely 
competent  or  there  is  more  to  the  production  of  software  for  the  mass  market 
than  meets  the  eye.  If  the  latter,  then  what  precisely  are  the  difficulties?  We  all 
await  a  succinct  description  of  the  intellectual  problems  that  would  establish 
mathematical  software  as  a  genuine  subdiscipline  of  computer  science. 

hast,  but  not  least,  we  should  emphasize  the  effort  expended  by  the  team 
on  the  documentation  and  uniform  coding  style.  It  is  a  mistake  to  speak  of 
"good"  documentation  because  documentation  which  satisfies  A  may  not  be 
suitable  for  B  .  It  is  reasonable  to  speak  of  documentation  being  suitable  for  a 
given  class  of  users.  The  EISPACK  GUIDE  has  intimidated  some  users  but.  at  the 
time  it  appeared,  it  glowed  with  virtue  in  comparison  with  other  examples  of 
documentation.  Moreover  nothing  on  this  scale  bad  been  attempted  before  by 
the  experts  on  mathematical  software. 


Inadvertently  EISPACK  has  provided  a  common  vocabulary  (the  names  of  its 
subroutines)  to  eigenvalue  hunters.  EISPACK  was  a  highly  visible  distillation  of 
what  had  been  learnt  by  the  matrix  eigenvalue  community  during  the  previous 
20  years.  It  was  good  for  public  relations.  It  predated  UNPACK  by  3  or  4  years 
and  yet  most  people  consider  the  eigenvalue  problem  to  be  significantly  more 
difficult  than  the  linear  equations  problem. 

To  numerical  analysts  EISPACK  seemed  to  be  the  solution  to  the  practical 
eigenvalue  problem.  What  else  remained? 

1.  EISPACK  routines  are  not  fast  enough  for  so-called  real  time  computa¬ 
tions  where  the  output  is  needed  in  microseconds. 

2.  EISPACK  manipulates  matrices  as  conventional  two  dimensional  arrays 
or  uses  a  conventional  one  dimensional  representation  of  symmetric  matrices. 
Except  for  the  tridiagonal  and  Hessenberg  form.  EISPACK  does  not  try  to  exploit 
zero  elements  in  the  original  matrix.  Moreover  the  stable  similarity  transforma¬ 
tions  used  by  most  of  the  subroutines  will  quickly  destroy  any  initial  sparsity. 

So  the  field  is  not  dead  by  any  means.  As  users  become  more  sophisticated 
they  are  finding  more  and  more  need  for  eigenvalues  and  eigenvectors.  So  usage 
grows. 

As  the  fast  storage  capacity  of  many  computer  systems  continues  to  grow 
so  does  the  domain  of  EISPACK.  Sound  advice  to  a  casual  user  is  to  use  EISPACK 
whenever  possible,  even  if  the  matrix  is  500  by  500  and  sparse. 

Who  will  then  remain  unsatisfied?  The  next  section  supplies  an  answer. 

3.  WHO  WANTS  EIGENVALUES  OF  LARGE.  SPARSE  MATRICES? 

EISPACK  subroutines  have  been  used  quite  extensively  and  it  has  been 
assumed  that  the  population  of  users  is  so  large  and  so  diverse  that  there  is  lit¬ 
tle  point  In  examining  the  market  more  closely.  Nevertheless  it  would  be  nice  to 
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know  how  strong  the  demand  is  for  programs  which  compute  some  eigenvalues 
and  eigenvectors  of  a  small  nonsymmetric  matrix. 

Even  less  warranted  is  the  natural  assumption  that  there  is  a  general  need 
for  an  extension  of  EISPACK  to  sparse  -  and  large  -  matrices.  Conversations  with 
a  variety  of  users  over  several  years  have  forced  the  author  to  the  following 
surprising  explanation  of  the  current  state  of  affairs. 

The  market  for  eigenvalue  programs  can  be  divided  into  two  quite  different 
camps:  INTENSIVE  users  and  SPORADIC  users. 

The  intensive  users  already  spend  millions  of  dollars  a  year  on  the  extrac¬ 
tion  of  eigenvalues  because  spectral  analysis  is  essential  to  their  daily  work. 
They  need  efficiency  and  have  already  crafted  their  programs  to  exploit  the  spe¬ 
cial  features  of  their  tasks.  General  purpose  programs  with  meticulous  code 
designed  to  cope  with  any  difficulty  which  might  arise  are  not  cost  effective  for 
this  group.  Of  course,  the  more  enlightened  intensive  users  will  collaborate  with 
experts  as  their  special  purpose  software  evolves  to  meet  even  more  exacting 
demands.  Actually  the  class  of  intensive  eigenvalue  hunters  also  splits  into  two 
quite  distinct  subclasses:  those  with  small,  dense  matrices  who  need  the  output 
in  real  time  (microseconds)  and  those  who  generate  larger  and  larger  matrices 
as  they  make  their  mathematical  models  more  realistic.  The  "real  time”  users 
may  be  driven  to  solve  their  problems  wuh  hardware  rather  than  software  and 
we  will  concentrate  here  on  the  large  rr.r.r.x  problem. 

The  SPORADIC  user  has  neither  the  incentive  nor  the  inclination  to  study 
matrix  methods.  The  need  for  some  eigenvalues  arises  and  the  user  wants  to 
obtain  them  with  minimal  fuss.  Reliability  is  more  important  than  efficiency  and 
EISPACK  is  the  answer  to  his  prayers.  We  suggest  that  the  number  of  sporadic 
users  whose  matrices  are  too  large  for  EISPACK  is,  and  will  remain,  very  low. 


The  foregoing  remarks  do  not  lessen  the  value  of  developing  good  methods 
for  all  sorts  of  large  eigenvalue  problems.  On  the  other  hand  the  situation 
should  make  one  hesitate  before  undertaking  the  painful  chore  of  producing  a 
general  purpose  reliable  package  for  most  sparse  eigenvalue  problems. 

INTENSIVE  USERS 

The  author  would  welcome  corrections  and/or  extensions  to  the  following 

list. 

I.  Structural  Engineers 

Most  eigenvalue  calculations  arise  at  the  heart  of  analyses  of  structures 
subject  to  small  vibrations.  There  is  an  n  x  n  global  stiffness  matrix  K  and  an 
n  x  n  global  mass  matrix  U  .  Both  are  symmetric  and  real.  K  is  positive 
definite.  The  usual  task  is  to  find  all  the  eigenvalues  ^  in  a  given  interval  at 
the  left  end  of  the  spectrum  (i.e.  near  0).  together  with  their  eigenvectors  z*  , 

(K-\  M)zi  -  0,i  =  l,2,  •  •  • 

The  z{  determine  the  shapes  of  the  fundamental  modes  of  vibration  of  the 
structure  and  the  \  determine  the  associated  natural  frequencies  2n/  ,  i 

”  1,2,  ...  . 

The  nonzero  elements  of  K  and  JJ  are  integrals  and  more  arithmetic 
operations  are  needed  to  form  K  and  7  than  to  compute  a  few  eigenvalues. 

In  1978  one  company  paid  12.0C,;  dollars  to  obtain  30  eigenvalue/vector 
pairs  for  a  problem  with  n  =  12,000.  This  cost  excludes  program  development. 
A  good  finite  element  package  was  used  with  the  technique  called  subspace 
iteration  for  the  eigenvalue  extraction.  Professor  E.  Wilson  (Civil  Engineering 
Department,  University  of  California,  Berkeley)  estimates  that  structural 
engineers  spent  about  10  million  dollars  in  1978  on  eigenvalue  computations.  A 
typical  problem  today  will  have  n  =  500  and  20  modes  computed. 


Nevertheless,  there  is  continued  demand  to  analyze  more  complicated  struc¬ 
tures  in  greater  detail.  Problems  with  n  >  10,000  and  400  modes  have  been 
solved. 

D.  Quantum  Chemists 

The  method  of  Configuration  Interaction  (Cl)  has  become  the  preferred 
method  for  deducing  observable  properties  of  real  and  hypothetical  molecules 
from  their  detailed  electronic  stucture.  The  Cl  method  approximates  the  solu¬ 
tion  of  the  clamped  nuclei  Schroedinger  equation  by  expanding  it  in  terms  of 
orthogonal  functions  made  out  of  products  of  so-called  single  and  double  elec¬ 
tron  spin  orbitals.  More  details  are  given  in  [Shavitt,  1977]  and  [Davidson,  1982]. 
These  papers  show  that  interesting,  difficult,  special  eigenvalue  problems  are 
being  solved  quite  independently  of  the  numerical  analysis  community.  Great 
ingenuity  goes  into  the  calculation  of  the  real  symmetric  matrix  H  (//for  Hamil¬ 
tonian  ).  The  nonzero  elements  of  H  are  multiple  integrals  and  constitute  only 
10%  of  the  positions.  Unfortunately  they  are  scattered  over  the  matrix  in  a  way 
that  precludes  any  simple  structural  template.  Each  eigenvector  represents  a 
wave  function  and  its  eigenvalue  is  the  energy  of  the  associated  state. 

In  these  chemical  computation*  the  determination  of  H  requires  10  to 
100  times  more  work  then  the  extraction  of  the  eigenvector  belonging  to  the 
smallest  (i.e.  most  negative)  eigenvalue.  Perhaps  100,000  dollars  are  spent  per 
year  in  the  U.S.A.  on  the  actual  eigenvalue/vector  computation.  Usually  only 
the  smallest  pair  is  required.  A  typical  calculation  has  the  order  n  =  10,000 
but  this  activity  is  expected  to  increase  sharply.  During  the  summer  of  1982  the 


ethylene  molecule  C2//4  was  analyzed  in  joint  work  at  Cambridge  University 
and  the  University  of  California,  Berkeley.  For  this  problem 

n  =  1,046,758 

This  calculation  demonstrated  convincingly  the  need  to  include  excited  states  in 
the  expansion.  In  other  words  the  simpler  models  in  Cl  are  not  adequate  for  the 
detail  required  in  current  investigations. 

It  is  sad  that  most  numerical  analysts  have  never  heard  of  the  eigenvalue 
methods  invented  by  the  chemists.  These  are  described  in  [Davidson,  1962].  The 
general  techniques  favored  by  the  numerical  analysts  are  too  inefficient  for  seri¬ 
ous  (i.e.  specialized)  tasks. 

I  have  not  identified  a  third  group.  In  the  1960's  there  was  considerable 
effort  expended  on  the  calculation  of  the  criticality  constant  for  nuclear  reac¬ 
tors  but  that  activity  has  subsided.  There  is  much  interest  in  the  vibration  of 
the  piping  systems  in  modern  reactors  but  that  work  is  part  of  structural 
analysis.  One  candidate  for  the  third  position  is  circuit  analysis  and  control 
problems  but  at  present  that  group  seems  to  favor  direct  solution  of  their  non¬ 
linear  differential  equations. 

The  intensive  users  have  developed  their  own  eigenvalue  software.  In  fact 
the  structural  engineers  discovered  the  method  of  simultaneous  iteration  for 
themselves  but  gave  it  a  more  descriptive  name,  subspace  iteration.  The 
method  itself  is  quite  obvious  and  what  counts  is  the  implementation.  It  is  sad 
that  the  beautiful  program  RITZIT,  developed  by  H..  Rotishauser  in  1968/69  and 
published  in  the  handbook  of  Wilkinson  and  Reinsch,  had  no  influence  on  the 
structural  engineers  working  in  the  U.S.A.,  although  in  Britain  the  work  of  Jen¬ 
nings  did  employ  some  of  the  techniques  embodied  in  RITZIT.  By  the  time  the 
Rii'UT  quality  does  creep  into  the  Unite  element  packages  then  subspace  itera¬ 
tion  will  have  been  displaced  by  modern  versions  of  the  Lanczoe  algorithm.  The 
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story  is  sad  because  numerical  analysts  should  justify  their  professional  support 
by  ftnmmnnlimtim  appropriate  techniques  to  users.  It  is  not  sufficient  to 
develop,  analyze  and  publish.  Moreover,  for  reasons  both  good  and  bad,  the 
engineers  will  not  read  the  numerical  analysis  journals  and  so  the  missionary 
wcrk  will  have  to  be  done  by  word  of  mouth. 

All  the  knowledge  acquired  by  the  specialists  in  matrix  computations  played 
perhaps  no  role  at  all  in  90%  of  the  sparse,  large  eigenvalue  extractions  made  in 
1975  in  the  U.S.A.,  i.e.  in  the  calculations  of  the  structural  engineers  and  the 
theoretical  chemists.  By  1985  the  situation  may  have  been  reversed.  Numerical 
analysts  such  as  P.S.  Jensen  at  Lockheed  (Palo  Alto)  and  J.  Lewis  at  Boeing  Com¬ 
puter  Services  (Seattle)  and  their  colleagues  are  deeply  involved  in  ambitious 
structural  analyses  and  mil  adapt  the  best  techniques  they  can  find  to  produce 
really  effective  software.  Portability  is  nice,  it"  is  particularly  important  for 
packages  aimed  at  general  users,  but  let  us  hope  that  it  never  becomes  a  com¬ 
mandment.  The  glamorous  work,  at  the  leading  edge  of  our  capacities,  will  have 
to  exploit  every  feature  that  is  special  to  the  problem  and  the  computer  system. 
This  is  the  way  of  Mary  (see  the  next  paragraph  if  you  are  not  familiar  with  this 
expression). 

Martha  complained  to  Jesus  that  while  she  was  busy  preparing  the  meals 
her  sister  Mary  just  listened  to  him.  Jesus  condoned  Mary's  behavior  and,  in 
some  traditions,  the  way  of  Martha  has  stood  for  the  worthy,  essential,  but  dull 
chores  (editorial  work,  perhaps?).  In  contrast  the  way  of  Mary  is  "where  the 
action  is". 

Those  sporadic  users  for  whom  EISPACK  is  inadequate  and  who  cannot 
develop  programs  for  themselves,  seem  to  be  in  hiding.  We  presume  they  exist 
and  so  the  mundane  but  worthy  task  of  providing  robust,  easy-to-use,  portable 
software  for  standard  sparse  problem  types  devolves  on  the  matrix  specialists. 
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11  is  the  way  of  Martha.  We  hope  that  the  reader  will  forgive  allusion  to  the  New 
Testament  [Luke  10:38-42]. 

5.  HOW  TO  EXPLOIT  SPARSITY. 

It  is  possible  to  approximate  an  eigenvalue  of  a  linear  operator  by  "sam¬ 
pling"  at  well  chosen  points  in  its  domain.  Thus  the  simple  power  method  sam¬ 
ples  the  action  of  a  matrix  A  on  the  sequence  of  column  vectors 
x.  Ax,  A*x,  Aax,  ■  •  The  sequence  of  Rayleigh  quotients  of  these  vectors  con¬ 
verges  quite  rapidly  to  the  dominant  eigenvalue  of  A  .  Recall  that  the  Rayleigh 
quotient  of  a  vector  v  is  vtAu/vlv  .  Consequently  there  are  methods  that 
need  only  be  given  a  subroutine,  call  it  OP  ,  which  returns  Atu  for  any  given 
vector  v . 

This  is  the  only  way  in  which  A  appears  in  the  method.  The  structure  and 
sparsity  of  A  can  be  exploited  by  the  user  in  writing  code  for  OP .  In  other 
words,  the  buck  is  passed  to  the  user!  He,  or  she.  is  in  the  best  position  to  take 
advantage  of  As  characteristics  to  speed  up  the  formation  of  Au  .  When  v  has 
10,000  components  one  pays  attention  to  this  product.  All  the  software  we 
describe  below  actually  ignores  sparsity  completely. 

This  is  in  stark  contrast  to  the  direct  solution  of  linear  equations  where  a 
number  of  clever  devices  are  used  to  take  advantage  of  zero  elements.  Some  of 
those  sparse  techniques  may  be  used  by  the  subprogram  OP  ,  but  none  of  that 
work  appears  explicitly  in  the  eigenvalue  codes.  Nevertheless  the  development 
of  methods  based  on  the  user-supplied  subprogram  OP  represents  a  beautiful 
division  of  labor  the  method  is  not  obscured  by  details  of  As  structure. 

This  is  the  place  to  mention  a  serious  confusion  that  has  arisen  in  the  past 
in  the  assessment  of  methods.  We  focus  on  symmetric  matrices  for  the  rest  of 
this  section.  At  a  certain  level  of  abstraction  the  power  method  is  identical  to 
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inverse  iteration.  One  technique  works  with  A  ,  the  other  with  A-1  .  Inverse 
iteration  performs  the  useful  task  of  finding  the  eigenvalue  nearest  to  0  but 
pays  the  significant  price  of  factoring  A  (toLU)  in  order  to  form  u  =  A~lv  (by 
solving  Lai  =  v  and  Uo  -  u)  .  If  a  sparse  eigenvalue  program  is  used  to  com¬ 
pute  eigenvalues  of  A  it  makes  an  ENORMOUS  difference  whether  the  subroutine 
OP  delivers  Au  or  A-1 v  . 

Subspace  iteration  finds  eigenvalues  close  to  a  by  factoring  A  —  al  anc. 
letting  OP  deliver  (A  -  al)~lv  .  When  a  =  0  the  eigenvalues  near  0  are  com¬ 
puted  first  and  quite  rapidly. 

The  Lanczos  algorithm  is  somewhat  different  It  produces  eigenvalues  at. 
both  ends  of  the  spectrum  of  the  linear  operator  represented  by  OP .  the  more 
extreme  ones  emerging  before  the  interior  ones  .  Thus  Lanczos  offers  the  hope 
of  computing  the  eigenvalues  of  A  nearest  0  without  the  pain  of  invoking  some 
sparse  factorization  of  A  .  In  this  sense  Lanczos  has  been  compared  with  Sub¬ 
space  Iteration.  The  performance  depends  quite  strongly  on  the  matrix  but  as 
the  order  n  grows  (n  >  100  say)  Lanczos  fares  worse  and  worse  and  is  soon 
eclipsed  by  Subspace  Iteration.  Lanczos  will  have  computed  perhaps  50 
unwanted  eigenvalues  near  «°  for  every  eigenvalue  near  0  .  This  is  because,  in 
most  given  problems,  the  larger  eigenvalues  are  much  better  separated  than 
the  small  ones.  (OP  is  approximating  an  unbounded  operator.)  Although  Lanc¬ 
zos  is  optimal  in  certain  important  respects  it  is  a  hopeless  task  to  compute  the 
smallest  eigenvalue  without  either  factoring  A  ,  or  solving  Au  -  6  iteratively, 
or  having  a  good  approximation  to  the  wanted  eigenvector.  This  is  certainly  the 
case  when  n  >  1,000  . 

The  solution  is  easy.  Give  Lanczos  the  same  OP  subroutine  as  Subspace 
Iteration.  Then  the  power  of  Lanczos  reveals  itself  quickly.  Both  methods  are 
then  working  with  (A  -  a)"1  to  find  eigenvalues  near  e  . 


8.  SOFTWARE  TOR  THE  HASSES 


By  the  time  a  matrix  expert  has  seen  his  program  appear  in  a  referreed 
journal  he  probably  never  wants  to  see  it  again.  The  program  comes  to  dom¬ 
inate  its  creator.  Dr.  R.C.  Ward  (ORNL)  is  well  aware  of  the  consequent  difficulty 
of  getting  good  matrix  programs  into  the  hands  of  non-expert  users.  In  conjunc¬ 
tion  with  the  Sparse  Matrix  Symposium  of  1982  in  Fairfield  Glade.  Ward,  and  his 
colleagues  at  Oak  Ridge,  have  started  a  public  catalog  of  software  for  sparse 
matrix  problems,  including  eigenvalue  extraction.  This  will  help  to  focus  the 
production  of  good  software. 

A  non-expert  will  still  be  annoyed  at  seeing  perhaps  six  programs  all 
designed  for  the  same  task  and  all  based  on  say  the  Lanczos  algorithm.  His 
annoyance  is  forgiveable  but  unwarranted.  He  should  be  made  to  understand 
that  the  sparse  eigenvalue  problem  is  still  a  research  topic.  The  experts  do  not 
yet  know  the  best  way  to  implement  the  Lanczos  algorithm,  (to  block  or  not  to 
block,  to  orthogonalize  the  Lanczos  vectors  a  lot,  a  little,  or  not  at  all)  or  how  to 
handle  secondary  storage.  Thus  it  is  good  that  there  are  several  rival  programs 
until  a  reasonable  consensus  is  reached.  Too  many  programs  is  better  than 
none  at  alL  For  this  reason  there  is  no  shame  in  the  fact  that  the  software  avail¬ 
able  now  is  not  up  to  E1SPACK  standards.  In  particular  none  of  the  codes  has 
been  subjected  to  widespread  testing  in  a  variety  of  computer  systems.  It  is 
worth  mentioning  that  Professc-  Ruhe  had  some  difficulty  in  transferring  his 
code  STLM  from  the  CDC  6600  or.  which  it  was  developd  in  Sweden  to  an  identical 
machine  at  Boeing.  His  code  suddenly  became  much  less  efficient.  Why?  The 
operating  system  was  different!  The  good  ship  Portability  may  well  founder  on 
the  rock  called  Operating  System. 

All  the  programs  mentioned  below  compute  one  or  more  eigenvaiue/vector 
pairs  of  the  symmetric  problem  ( A  -  \B)*  =  0  ,  unless  the  contrary  is  stated. 


comments  are  given  at  the  side.  All  programs  are  in  FORTRAN,  most  in  the  ANSI 
standard  version  of  FORTRAN  66.  All  are  portable  to  some  extent,  some  are  com¬ 
pletely  portable.  Between  30  and  60%  of  the  lines  are  COMMENTS.  A  potential 
user  should  consult  Ward's  catalog  for  more  information  on  the  programs  men¬ 
tioned  below. 

A.  Programs  based  on  Lanczos 

The  modern  versions  of  the  Lanczos  algorithm  are  not  simple  and  most  of 
them  are  described  in  Chapter  13  of  [Parlett.  1980].  A  brief  summary  will  be 
given  here,  but  please  see  the  discussion  in  Section  8  for  some  important 
details. 

Lanczos  algorithms  are  iterative  in  nature.  A  starting  vector  (or  block  of 
vectors)  is  chosen  and  at  each  step  a  new  vector  is  added  to  the  sequence.  In 
exact  arithmetic  these  vectors  are  mutually  orthogonal  and  of  unit  length.  A 
characteristic  and  pleasing  feature  is  that  only  the  two  most  recent  vectors  are 
needed  in  the  computation  of  the  next  one.  In  addition  to  these  vectors  the 
Lanczos  method  builds  up  a  symmetric  tridiagonal  matrix  T .  each  step  adds  a 
new  row  (and  column)  to  T.  Quite  soon  some  eigenvalues  of  T  begin  to  approx¬ 
imate  some  eigenvalues  of  the  operator  hidden  in  OP .  Almost  always  it  is  the 
extreme  eigenvalues  which  are  well  approximated. 

Sometimes  an  extreme  eigenvalue  is  approximated  to  full  working  precision 
(say  15  decimals)  after  only  30  Lanczos  steps.  It  is  rare  that  one  can  solve  a 
linear  system  correct  to  4  deciamls  after  only  30  stpes  of  the  conjugate  gra¬ 
dient  algorithm  (which  is  intimately  related  to  the  Lanczos  algrithm),  unless  an 
excellent  preconditioner  is  used.  Thus  it  seems  "easier"  to  And  a  few  extreme 
eigenvalues  and  eigenvectors  of  A  than  to  solve  Ax  *  b  \ 
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Tha  codas 

All  of  the  programs  In  Table  A  except  for  EA14,  stand  alone,  they  do  not 
invoice  other  packages  explicitly.  This  feature  turns  out  to  be  attractive  to  new 
customers  although  it  is  irritating  that  EISPACK  subroutines  must  be  copied  into 
the  program  instead  of  being  invoked. 

The  only  program  which  can  be  obtained,  independently  of  the  author,  from 
the  Argonne  Code  Center  (now  called  NESC,  the  National  Energy  Software 
Center)  is  Scott’s  LAS02.  Documentation  is  available  on  an  accompanying  file. 
This  code  has  been  used  extensively  at  Oak  Ridge  National  Laboratory  on  a 
variety  of  structural  analysis  problems. 

Documentation  for  the  Cullum  and  Willoughby  codes  is  obtainable  in  book 
form  (like  the  EISPACK  Guide).  Their  programs  use  little  or  no  reorthogonaliza- 
tion  of  the  Lanczos  vectors  but  have  developed  ingenious  ways  to  identify  which 
eigenvalues  of  the  auxiliary  tridiagonal  matrix  actually  belong  to  the  original 
matrix.  Their  code  is  shorter  than  the  rival  codes. 

The  Swedish  program  uses  blocks  but  does  little  reorthogonalization.  Its 
chief  feature  is  a  sophisticated  mechanism  for  choosing  the  sifts  a  in  the  spec¬ 
tral  transform  technique  launched  by  Ruhe  and  described  in  Section  8.  A  careful 
comparison  of  LAS02  with  STLM  would  be  very  interesting  for  the  small  band  of 
specialists  in  matrix  eigenvalue  computations  and  would  help  in  progress 
cowards  a  preferred  implementation  of  the  Lanczos  algorithm. 

The  program  EA14  is  much  shorter  than  the  others  because  it  does  not  use 
blocks,  does  not  do  any  reorthogonalization,  and  finds  eigenvalues  only,  and  not 
even  their  multiplicity.  The  user  selects  the  interval  to  be  explored.  If  the 
interval  happens  to  be  empty  the  code  will  report  that  fact  in  reasonable  time. 
This  is  noteworthy  because  the  code  assumes  that  OP  delivers  An  and  so  tri¬ 
angular  factorization  is  not  available.  Thus  the  standard  technique  for  checking 


the  number  of  eigenvalues  in  an  interval  is  not  available.  Ian  Gladwell  is  incor¬ 
porating  EA14  into  a  program  which  also  computes  eigenvectors.  EA15.  This  will 
also  go  into  the  NAG  library. 
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Table  A 

Symmetric  Problems,  Lanczos  Based. 


Author 

Name 

Lines 

Distributor 

Comments 

J.  Cullum  It 

LMAJN 

2100 

authors 

A  -  A/  only. 

R.  Willoughby 

LVMAIN 

(IBM) 

no  orthogonalization. 

J.  Cullum  It 

R  Willoughby 

BLMAIN 

1000 

H 

A  -  \I  only. 

block  limited  reorthog. 

D.S.  Scott 

LAS02 

3288 

NESC 

block. 

selective  orthog. 

T.  Ericsson 

It  A.  Ruhe 

STLM 

8500 

authors 

(Sweden) 

Shift  and  Invent 

strategy.  (See  Section ) 

no  orthogonalization. 

B.N.  Parlett 

It  J.  Reid 

EA14 

648 

Harwell 

All  eigenvalues  in 

a  given  interval. 

No  eigenvectors. 

No  orthogonalization. 

A  -XI  only. 
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B.  Programs  not  based  on  Lancsos. 

There  are  a  number  of  long,  sophisticated,  implementations  of  subspace 
iteration  buried  in  Finite  Element  packages.  As  such  they  are  beyond  the  limits 
of  survey.  However  the  first  three  programs  in  Table  B  are  available  realiza¬ 
tions  of  subspace  iteration. 

This  paragraph  sketches  the  method.  An  initial  set  of  orthogonal  vectors  is 
chosen.  The  number  of  vectors  in  the  set  is  the  dimension  of  the  subspace.  Call 
it  p  .  The  proper  choice  of  p  is  impor  tant  and  difficult.  These  vectors  may  be 
considered  as  the  columns  of  an  nxp  matrix  X° .  At  step  j  the  subroutine  OP 
is  used  to  compute  Y  :=  (.4  -  oB)~l  X*-1  .  Next  the  Rayleigb-Ritz  approxima¬ 
tions  from  the  column  space  at  Y  are  computed.  These  Ritz  vectors  go  into  the 
columns  of  the  matrix  Xi  .  They  provide  the  best  orthonormal  set  of  approxi¬ 
mate  eigenvectors  that  can  be  obtained  by  taking  linear  combinations  of  the 
columns  of  Y .  After  a  while  some  of  the  columns  of  X*  provide  excellent 
approximations  to  some  eigenvectors  of  the  pair  (A.P)  .  There  are  a  number  of 
variants  of  this  scheme  and  several  clever  tricks  in  the  implementation.  See 
[Parlett,  1960,  Chap.  14],  [Jennings.  1981]  or  [Stewart,  1976]  for  more  details. 

The  first  three  programs  in  Table  B  are  realizations  of  subspace  iteration. 
The  Achilles  heel  of  the  technique  is  the  selection  of  block  size.  The  first  two 
programs  are  by  numerical  analysts  and  it  would  be  intersting  to  see  how  they 
compare  with  the  codes  imbedded  in  the  finite  element  (FEM)  packages  men¬ 
tioned  in  Section  2.  It  is  likely  that  they  are  shorter,  cleaner,  more  robust,  and 
more  efficient  than  their  FEM  rivals.  On  the  other  hand  they  have  not  been  fine 
tuned  for  structural  analysis  problems  and  that  might  make  a  difference. 

The  entry  TRACMN  is  the  product  of  recent  research.  It  is  based  on  the  fact 

frocs  [F*(4  **'  oB)F\ 


that 
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is  minimized  over  all  n  x  n  matrices  F  when,  and  only  when,  the  columns  of 
F  are  the  eigenvectors  of  F  belonging  to  the  p  eigenvalues  closest  to  a.  At 
each  step,  the  current  F  is  adjusted  in  order  to  reduce  the  trace.  We  hope  that 
the  authors  will  compare  it  with  LAS02  or  some  similar  rival  technique. 

L0PS1  is  the  only  program  offered  in  Ward's  catalog  for  the  nonsymmetric 
problem.  The  program  has  been  published  in  TOMS  (a  severe  test)  and  employs 
subspace  iteration. 

One  version  of  Davidson's  method  is  realized  in  DIAG  and  was  listed  in  the 
catalog  of  the  now  defunct  National  Resource  for  Computation  in  Chemistry 
(LBL,  Berkeley).  The  general  reader  is  warned  that  Davidson's  method  has  been 
developed  for  problems  in  Quantum  Chemistry.  The  matrix  must  be  strongly 
diagonally  dominant  in  order  for  its  perturbation  technique  to  be  justified.  When 
all  the  diagonal  elements  of  H  are  the  same  then  Davidson's  method  reduces  to 
Lanczos.  The  difference  is  as  follows.  At  step  j  the  Rayleigh-Ritz  technique  is 
used  to  produce  the  best:  approximation  Vj  to  the  fundamental  eigenvector 
that  can  be  obtained  as  a  linear  combination  of  the  current  basis  vectors 
6lt  ....  bj  .  Let  p(x)  be  the  Rayleigh  quotient,  p(z)  =  xlHx/zlx  .  The  residual 
vector 

rj  ■  -  (H  -p(l/j))Vy 

is  proportional  to  the  gradient  vector  Vp  at  i;  .  > 


Perturbation  theory  now  says  that  if  v;  is  a  good  approximation  then  an 
even  better  one  is 

*j  '•  =  (p(vj)  ~  diag(H))-1  rj 

Davidson  proposed  to  take  as  bj+ 1  the  normalized  component  of  ttj  orthogonal 
to  6  j . bj  . 


We  make  three  observations.  1.  The  "interaction"  matrix  or  projection 
B*  HB  is  not  tridiagonai  but  full.  2.  The  method  aims  at  a  particular 
eigenvalue/vector  pair.  3.  When  diag(H)  =  hnl  then  is  a  multiple  of  r}  and 
Lanczos  is  recovered.  This  last  assertion  is  not  obvious,  since  Vj  *  bj  ,  but  it  is 
true. 

Observation  2  shows  that  Davidson  is  not  really  a  rival  to  Lanczos.  The 
methocs  address  different  tasks.  Moreover  Davidson  exploits  the  face  that  good 
starting  vectors  are  available  from  chemistry. 


Table  B 

Symmetric  Problems,  other  methods. 

None  of  the  programs  stands  alone  (TOMS  =  published  in 
the  Transactions  on  Mathematical  Software  of  the  ACM). 


Author 

Name 

Lines 

Distributor 

Comments 

1.  Duff 

EA12 

427 

Harwell 

(England) 

A  -  A/  only. 

RITZIT  inspired 

P.J.  Nikolai 

SIMITZ 

550 

IMSL 

author 

A  —\B 

RITZIT  inspired 
(in  TOMS) 

A.  Jennings 

SI 

7 

author 

(Ireland) 

A  —\B 

(in  Int.  Joum.  Num. 
Methods  in  Engineering) 

J.A.  Wisniewski 
A.H.  Sameh 

TRACMN 

5B6 

authors 

A  —  \B 
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7.  LANCZOS  VERSUS  SUBSPACE  ITERATION 

These  two  rival  techniques  are  good  and  address  the  same  task.  A  com¬ 
parison  of  them  is  given  in  [Parlett,  Nour-Omid,  Taylor  1983]  and  the  conclusion 
there  is  that  a  modern,  properly  used  version  of  the  Lanczos  algorithm  should 
be  an  order  of  magnitude  (base  10)  faster  than  a  good  subspace  iteration  pro¬ 
gram  on  problems  with  n  >  500  .  The  larger  the  number  of  modes  wanted  the 
greater  the  advantage  of  the  Lanczos  algorithm. 

The  fundamental  theoretical  advantage  of  the  Lanczos  is  that  it  never  dis¬ 
cards  any  information  whereas  subspace  iteration,  at  each  step,  overwrites  a  set 
of  approximate  eigenvectors  with  a  better  set.  What  is  surprising  is  that  the 
simple  Lanczos  algorithm  needs  only  5  or  6  n-vectors  in  the  fast  store  at  each 
step.  Moreover,  it  can  find  multiple  eigenvalues  if  selective  orthogonalization  is 
used. 

There  are  enough  poor  implementations  of  Lanczos  available  to  complicate 
comparisons.  The  engineers  who  have  devoted,  themselves  to  steady  improve¬ 
ments  of  subspace  iteration  are  loyally  defending  the  virtues  of  their  codes.  It 
will  be  interesting  to  see  what  happens. 

8.  WHAT  HAVE  WE  LEARNT  SINCE  1978? 

Much  could  be  said  in  answer  to  this  question  but  two  items  suggest  them¬ 
selves. 

1.  The  importance  to  new  users  of  stand  alone  programs.  This  point  was 
made  in  Section  5. 

2.  The  Spectral  Transformation.  The  Lanczos  algorithm  requires  that  a 
general  linear  problem 

[K  -  AAf  ]*  =  0 

be  reduced  to  standard  form.  There  are  several  ways  to  accomplish  this  and  it 


is  vital,  for  large  problems,  to  choose  a  good  one.  Although  the  facts  given  below 
are  elementary  and  have  been  known  to  the  specialists  in  matrix  computations 
for  a  long  time  it  is  fair  to  say  that  until  A.  Ruhe  pointed  out  their  implications  in 
[Ericsson  and  Ruhe  1980]  they  were  not  given  the  proper  emphasis.  Of  course, 
subspace  iteration  has  employed  spectral  transformations  from  the  earliest 
days. 

Ruhe  proposed  that  the  original  problem  be  rewritten  in  the  standard  form 
[MKK  -  <jlU)-'  M*  -  vl]y  ~  0 

with  y  -  M*x  ,  and  a  a  shift  into  the  interior  of  the  interval  which  is  to  be 
searched  for  eigenvalues  A  .  This  is  quite  different  from  the  usual  recommen¬ 
dation 

[L~'KL-1  -  A I]z  =  0 
where  M  -  LL*  and  *  =  L*z  . 

In  the  majority  of  large  structural  problems  (with  displacement  formula¬ 
tion)  M  is  singular  as  well  as  diagonal,  perhaps  1/3  of  its  diagonal  elements  are 
zero.  The  shift  a  is  not  fixed  precisely  so  that  there  is  no  loss  in  assuming  that 
K  —  <jM  is  nonsingular  and  can  be  factored  in  a  stable  manner  without  any  row 
and  column  interchanges.  If  K  -  oM  -  LDL *  then  the  product 
y  =  U #  ( K  -  oM)~l  M #  u  is  formed  in  stages 

v  *-  xl  ,  solve  L  w  =  v  .  sc. vs  DLlx  =  w  ,  y  «-  x  . 

If  the  factorization  of  K  -  oM  is  too  costly  then  in  principle  (AT  -  crM) y  =  M%u 
could  be  solved  by  an  iterative  method.  The  point  is  that  the  subprogram  OP  in 
Lanczos  should  return  M\K  -  u  when  given  tt  .  Sparse  techniques 

can  be  used  in  factoring  K  -  oM  .  This  reduction  of  the  problem  transforms  the 
spectrum  according  to 
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1 


Xt  -  a 

and  so  the  eigenvalues  \  closest  to  a  turn  into  the  eigenvalues  t/<  closest  to 
<>»  .  These  are  the  ones  which  will  appear  first  in  a  run  of  Lanczos  algorithm. 
For  large  n  the  reduction  in  the  number  of  steps  is  so  satisfactory  that  it  will 
offset  the  cost  of  factorization.  Without  the  spectral  transformation  it  will  often 
be  necessary  to  take  %  n  or  more  steps  of  the  Lanczos  algorithm  in  order  to 
approximate  the  smallest  10  eigenvalues.  Yet  the  Lanczos  algorithm  is  most 
effective  when  used  with  fewer  than  200  steps. 

9  WORK  IN  PROGRESS 

It  is  surprising,  at  first,  that  the  software  reviewed  above  is  so  narrowly 
focussed.  So  it  is  natural  that  a  few  investigators  are  working  to  "fill  the  gaps", 
to  provide  robust  programs  for  all  the  requests  that  might  conceivably  be  made 
for  eigenvalues  of  various  types  of  large  matrix. 

Normal  matrices  enjoy  a  full  set  of  orthonormal  eigenvectors  and  most 
methods  designed  for  real  symmetric  matrices  extend  naturally  to  normal  ones. 
The  most  important  subclass  of  normal  but  assymmetric  matrices  are  those 
that  are  skew  and  Ward  has  developed  efficient  software  for  them,  as  described 
in  [Ward.  1978]. 

Work  on  nonnormal  sparse  matrices  is  still  in  the  experimental  stage.  The 
experts  are  exploring  the  problem  and  there  is  no  concensus  on  a  preferred 
technique.  Y.  Saad  has  adapted  Arnoldi's  method  for  sparse  problems.  This  is 
the  natural  extension  of  the  Lanczos  idea  but  subject  to  the  constraint  of  using 
orthonormal  bases  for  the  associated  Krylov  subspaces.  The  auxiliary  matrix 
generated  in  the  course  of  the  algorithm  is  Hessenberg  (i.e  hq  =  0  if  i  >  j  +  l) 
rather  than  triadiagonaL  This  is  not  a  serious  feature  provided  that  fairly  short 
runs  of  Lanczos  are  used,  say  50  steps  at  most  Saad  is  experimenting  with 
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variations  on  this  method  in  which  orthogonality  is  given  up  in  return  for  a 
banded  structure  in  the  Hessenberg  matrix,  see  [Saad,  I960]. 

The  Lanczos  algorithm  was  generalized  to  the  nonnormal  case  by  Lanczos 
himself,  the  procedure  is  obvious.  Unfortunately  it  can  breakdown  and,  in  fact, 
does  so  more  often  than  not.  In  [Parlett  and  Taylor,  1981]  it  is  shown  how  to 
restore  a  large  measure  of  stability  in  return  for  modest  increase  in  complexity. 
Some  feature  of  the  original,  strict  Lanczos  algorithm  must  be  discarded  if  the 
breakdown  is  to  be  avoided.  The  new  method  looks  ahead  at  every  step  and 
decides  whether  to  enlarge  the  Krylov  subspace  by  one  or  two  dimensions.  The 
auxiliary  matrix  is  not  quite  triangular,  there  is  a  bulge  every  time  the  dimen¬ 
sion  increases  by  two.  Column  and  row  eigenvectors  are  given  equal  status  and 
condition  numbers  are  automatically  computed. 

Axel  Ruhe  is  experimenting  with  a  two-sided  Arnoldi  process  but  this  work  is 
at  an  early  stage. 

What  impedes  the  production  of  pleasing  software  for  nonnormal  problems 
is  the  potential  ill-condition  of  the  problem  itself.  Expectations  have  been  set  by 
the  symmetic  case  were  the  theory  is  most  satisfactory.  In  the  general  case  it  is 
difficult  to  generate  useful  error  estimates,  let  alone  error  bounds,  and  this 
affects  the  selection  of  stopping  criterion.  Numerical  analysts  would  want  to 
deliver  condition  numbers  along  with  the  eigenvalues,  but  users  do  not  want  this 
information,  especially  if  it  increases  the  cost  by  a  noticeable  amount. 

We  terminate  this  digression  with  a  reminder  that  Jenning's  program  LOPSI 
has  been  published  and  is  available  to  the  public.  LOPS!  adapts  subspace  itera¬ 
tion  to  the  asymmetric  case. 

The  most  promising  developments  are  close  to  the  symmetric  case.  Many 
problems  in  dynamics  arise  in  the  form 


where  x  is  the  time  derivative  of  the  vector  field  x{t )  .  The  damping  matrix  C 
is  such  a  pain  that  it  is  usually  ignored.  The  program  [Scott  St  Ward.  1982]  looks 
to  the  future  and  presents  software  for  the  associated  quadratic  eigenvalue 
problem 

(X2U  +  \C  +  K)u  =  0  . 

Their  method  is  based  on  the  Rayleigh  quotient  iteration.  It  does  not  take  the 
easy  way  out  and  reduce  the  quadratic  problem  to  a  linear  one  of  twice  the  size. 

The  Lanczos  programs  in  Ward’s  catalog  work  very  well  in  combination  with 
the  spectral  transformation  d  scussed  in  Section  8.  That  approach  is  not  possi¬ 
ble  when  the  user  cannot,  or  mil  not,  allow  the  solution  of  linear  systems  of  the 
form 

(K  -  ofi)v  =  u  . 

What  can  be  done? 

Scott  observed,  in  [Scott,  1981],  that  when  a  is  close  to  an  eigenvalue 
a  +  <J  of  the  pair  (KM)  them  the  vector  x  belonging  to  the  eigenvalue  5 
closest  to  zero  for  the  standard  problem 

(K  -a  U  -  kl)x  =  0 

is  an  approximate  eigenvector  of  (K.M)  belonging  to  a  .  He  formulates  an 
iterative  method  in  which  a  standard  sparse  eigenvalue  problem  is  solved  at 
each  step  for  (6.x)  and  a  +  <J  converges  monotonically  to  an  eigenvalue  of 
(KM)  . 


ICyONCLUSION 

The  eigenvalue  problem  for  large,  sparse  matrices  is  not  yet  understood 


well  eough,  in  all  its  computational  ramifications,  to  permit  the  rapid  deploy¬ 
ment  of  impeccable  software  to  the  masses  in  the  style  set  by  EISPACK.  Indeed 


it  is  doubtfull  that  the  public  is  impatient  for  the  arrival  of  this  facility.  Those 
with  an  urgent  need  for  such  software  prize  efficiency  above  generality.  They 
have  developed  their  own  programs  independently  of  the  numerical  analysts  and 
will  continue  to  do  so  unless  the  matrix  experts  go  to  them  and  demonstrate 
that  they  can  be  useful. 

The  wide  variety  of  computing  environments  is  going  to  play  havoc  with 
naieve  concepts  of  portabliy.  "Transmission  of  thought  is  one  of  the  two  basic 
themes  of  science.  The  interesting  and  difficult  task  is  to  find  the  appropriate 
level  at  which  to  transmit. 


******* 
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